Notes for cue vs uptake vs half-saturation trade-off invasion model
K Todd-Brown (ktoddbrown@gmail.com)
2017-02-14
Set up
library(assertthat)
library(reshape2)
library(ggplot2)
library(GGally)
library(plyr)
library(rootSolve)
library(deSolve)
library(knitr)
library(pander)
library(cowplot)
#sourceFiles <- 'R/invasion.R'
#l_ply(sourceFiles, source)
set.seed(100) #reproducible random selection
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X Yosemite 10.10.5
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cowplot_0.7.0 pander_0.6.0 knitr_1.15.1 deSolve_1.14
## [5] rootSolve_1.7 plyr_1.8.4 GGally_1.3.0 ggplot2_2.2.0
## [9] reshape2_1.4.2 assertthat_0.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.8 magrittr_1.5 munsell_0.4.3
## [4] colorspace_1.3-1 stringr_1.1.0 tools_3.3.2
## [7] grid_3.3.2 gtable_0.2.0 htmltools_0.3.5
## [10] yaml_2.1.14 lazyeval_0.2.0 rprojroot_1.1
## [13] digest_0.6.10 tibble_1.2 RColorBrewer_1.1-2
## [16] evaluate_0.10 rmarkdown_1.2 stringi_1.1.2
## [19] scales_0.4.1 backports_1.0.4 reshape_0.8.6Basic model description
The goal of this model to simulation biologically mediated decomposition. There are two main pools, Substrate (\(C\)) and Biomass (\(B\)). Carbon is added to the substrate pool via inputs (\(I\)), leached from both pools at a rate proportional to the carbon in the pool (\(h\) fraction of the biomass pool \(B\), and \(m\) fraction of the substate pool \(C\)), and transfer between the substrate \(C\) pool and biomass \(B\) pool via a Monod uptake with a maximum rate \(v_{max}\), half-saturation constant \(k\), and carbon use effiency \(\epsilon_{cue}\).
\(\frac{dB}{dt}=\frac{\epsilon_{cue}vBC}{k+C}-hB\)
\(\frac{dC}{dt}=I-mC-\frac{vBC}{k+C}\)
At steady state (\(\frac{dB}{dt}=\frac{dC}{dt}=0\)) this implies that:
\(B=\frac{\epsilon h}{I-mC}\)
\(C=\frac{hk}{\epsilon v-h}\)
Calculate the leaching and turnover rate
Assume that the input rate, biomass, and carbon pools are given targets, and that we are exploring cue, v, and k. We can thus define the leaching and turnover rates at steady state:
find_h <- function(C, cue, v, k){
return(C*cue*v/(k+C))
}
find_m <- function(cue, v, h, k, I, B){
return((cue*v-h)/(cue*k)*(cue/h*I-B))
}
If we further assume that there is a trade off between both 1) uptake \(v\) and carbon use effiency (\(\epsilon\)) such that \(v=\frac{e^{b\epsilon_cue}-e^b}{1-e^b}\) and the half saturation constant (\(k\)) such that \(k=a v +k_{min} = \frac{a}{1-e^b}(e^{b\epsilon_cue}-e^b)\). Defined as:
Visualize trade-off between cue vs uptake and k
cue_v_tradeoff <- function(b, vmax, cue){
return(vmax*(exp(b*cue)-exp(b))/(exp(0)-exp(b)))
}
tradeoff.df <- adply(.data=c(-0.1*1:9, -1*1:10, 0.1*1:9, 1:10), .margins=c(1), .id=c('id'), .fun=function(b){
ans <- data.frame(b=b, cue=seq(0, 1, length=1000))
ans$v <- cue_v_tradeoff(b=b, vmax=2, cue=ans$cue)
return(ans)
})
ggplot(tradeoff.df) + geom_line(aes(x=cue, y=v, group=b, color=b))
v_k_tradeoff <- function(kmin, k_v_slope, v){
return(kmin+k_v_slope*v)
}
tradeoff.df$k <- v_k_tradeoff(kmin=0, k_v_slope=100/2, tradeoff.df$v) #slope is max SOC over vmax rangeThis leads to the following steady state calculations for one biomass pool:
steadyState <- function(b, vmax, kmin, k_v_slope, h, m, I, cue){
v <- cue_v_tradeoff(b=b, vmax=vmax, cue)
k <- v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=v)
C <- h*k/(cue*v-h)
B <- cue/h*(I-m*C)
return(list(B=B, C=C))
}Steady state sensitivity to inputs
How do inputs affect stead state?
ssParm.ls <- list(b=c(-10, -1, 1, 10),
vmax=c(0.1, 1, 10, 100),
kmin=c(0, 1, 10, 100, 1000),
k_v_slope=c(1/10, 1, 10, 100, 1000),
h=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
m=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
cue=c(0.2, 0.4, 0.6, 0.8),
I=c(0.1, 1, 10))
steadyState.df <- expand.grid(ssParm.ls)
#Calculate the steady state values for each parameter combination
test <- with(steadyState.df, as.data.frame(steadyState(b, vmax, kmin, k_v_slope,
h, m, I, cue)))
steadyState.df <- cbind(steadyState.df, test)
#Check to see if it's a reasonable value over any inputs
steadyState.df$OnTarget <- with(steadyState.df, B+C > 10 & B+C < 500 &
B/(B+C)> 0.001 & B/(B+C) < 0.15)
goodPars <- ddply(steadyState.df, c('b', 'vmax', 'kmin', 'k_v_slope', 'h', 'm', 'cue'),
summarize, validSet=any(OnTarget))
goodPars <- subset(goodPars, validSet, -validSet)
goodPars$index <- 1:nrow(goodPars)
#Calculate the change in the steady state over different input values
ssGrad.df <- ddply(goodPars, c('b', 'vmax', 'kmin', 'k_v_slope', 'h', 'm', 'cue'),
function(xx){
xx <- data.frame(xx, I=c(0.1, 1, 10)*c(rep(1, 3), rnorm(3, mean=1, sd=0.05)))
temp <- cbind(xx, as.data.frame(with(xx, steadyState(b, vmax, kmin, k_v_slope, h, m,
I, cue))))
ans <- temp[1:3,]
ans$dB.dI <- (temp$B[1:3]-temp$B[4:6])/(temp$I[1:3]-temp$I[4:6])
ans$dC.dI <- (temp$C[1:3]-temp$C[4:6])/(temp$I[1:3]-temp$I[4:6])
return(ans)
})
ssGrad.df$OnTarget <- with(ssGrad.df, B+C > 10 & B+C < 500 & B/(B+C)> 0.001 & B/(B+C) < 0.15)
plot.df <- melt(subset(ssGrad.df, index %in% sample(unique(ssGrad.df$index), size=10),
select=c('index', 'I', 'B', 'C')),
id.vars=c('I', 'index'))
ggplot(plot.df) + geom_line(aes(x=I, y=value, color=as.factor(index))) +
facet_wrap(~variable) +
labs(title='Steady state carbon pool by inputs', x='Inputs', y='Pool size') +
guides(color=FALSE)
plot.df <- melt(subset(ssGrad.df, OnTarget),
id.vars=c('B', 'C', 'dB.dI', 'dC.dI', 'OnTarget', 'index'))
ggplot(plot.df) + geom_boxplot(aes(x=as.factor(value), y=dB.dI)) +
facet_wrap(~variable, scales='free') +
labs(title='dB/dI', x='Parameter value', y='Gradient (dB/dI)')
ggplot(plot.df) + geom_boxplot(aes(x=as.factor(value), y=dC.dI)) +
facet_wrap(~variable, scales='free') +
labs(title='dC/dI', x='Parmaeter value', y='Gradient (dC/dI)')
Optimum carbon use effiency
We want the optimum carbon use effiency where biomass is maximum as steady state but that is not possible to solve analytically (you end up with \(y=xe^x\), solve for x).
optimum_cue <- function(b, vmax, kmin, k_v_slope, h, m, I){
cue <- seq(0, 1, length=100)
ans <- steadyState(b, vmax, kmin, k_v_slope, h, m, I, cue)
ans$B[ans$B <= 0 | ans$C <= 0] <- -Inf
best <- which.max(ans$B)
return(list(B=ans$B[best], C=ans$C[best], cue=cue[best]))
}Generate reasonable parameter sets without compitition
Given that we have input targets between 0.1 and 10 mg-C per g-soil per day, total carbon pools of between 10 to 500 mg-C per g-soil, and a biomass to total carbon ratio of between 0.1 and 15 percent. Let’s generate some reasonable parameters to draw on.
parm.ls <- list(b=c(-10, -1, -0.1, 0.1, 1, 10),
vmax=c(0.1, 0.5, 1, 5, 10, 100),
kmin=c(0, 1, 10, 100, 1000),
k_v_slope=c(1/10, 1, 10, 100, 1000),
h=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
m=c(1e-4, 1e-3, 1e-2, 0.1, 1, 10),
I=c(0.1, 1, 10))
optParm.df <- expand.grid(parm.ls)
optParm.df <- ddply(optParm.df, names(optParm.df), function(xx){
ans <- as.data.frame(optimum_cue(b=xx$b, vmax=xx$vmax, kmin=xx$kmin, k_v_slope=xx$k_v_slope,
h=xx$h, m=xx$m, I=xx$I))
return(ans)
})
optParm.df$OnTarget <- with(optParm.df, B+C > 10 & B+C < 500 & B/(B+C)> 0.001 & B/(B+C) < 0.15)
#parm.allI <- ddply(subset(optParm.df, OnTarget), setdiff(names(parm.ls), 'I'), summarize,
# numI=length(I), minI=min(I), maxI=max(I))
#ggpairs(subset(optParm.df[,c(names(parm.ls), 'OnTarget')])
ggplot(subset(optParm.df, OnTarget)) + geom_histogram(aes(x=cue))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(subset(optParm.df, OnTarget)) + geom_violin(aes(x=as.factor(I), y=cue)) 
#ggplot(subset(optParm.df, OnTarget)) + geom_jitter(aes(x=I, y=B+C)) + scale_x_log10()Competition model co-existance
Consider the senerio where the native popuation is invated by a second competiting population that has a different carbon use effiency but is otherwise equivalent. Let the strategic carbon use effiency be where the native popuation is never invadable no matter what carbon use effiency value the invader has.
In otherwords both popuations are stable where one is either 0 or \(\frac{\epsilon_n v_n}{k_n + C} = \frac{\epsilon_i v_i}{k_i + C}\) where \(n\) and \(i\) stand for the native and invading populations respectively.
Again, not solvable analytically given our trade-off function but numerically we can calculate the following:
strategic_cueI <- function(b, vmax, kmin, k_v_slope, m, h, I){
cue <- seq(0, 1, length=100)
v <- cue_v_tradeoff(b=b, vmax=vmax, cue)
k <- v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=v)
ans <- steadyState(b, vmax, kmin, k_v_slope, h, m, I, cue)
measure <- cue*v/(k+ans$C)
measure[measure <= 0] <- NA
flag <- which.max(measure)
return(list(cue_strat=cue[flag], C_strat=ans$C[flag], B_strat=ans$B[flag]))
}
stratParm.df <- optParm.df
stratParm.df <- rename(stratParm.df, c('cue'='cue_opt', 'B'='B_opt', 'C'='C_opt', 'OnTarget'='validPools_opt'))
stratParm.df <- ddply(stratParm.df, names(stratParm.df), function(xx){
return(as.data.frame(strategic_cueI(b=xx$b, vmax=xx$vmax, kmin=xx$kmin,
k_v_slope=xx$k_v_slope, m=xx$m, h=xx$h, I=xx$I)))
})
stratParm.df$validPools_strat <- with(stratParm.df, B_strat+C_strat > 10 &
B_strat+C_strat < 500 &
B_strat/(B_strat+C_strat)> 0.001 &
B_strat/(B_strat+C_strat) < 0.15)
stratParm.df$orderDiff_mh <- round(log(stratParm.df$m)-log(stratParm.df$h), 1)
ggplot(subset(stratParm.df, validPools_opt & validPools_strat)) +
geom_point(aes(x=cue_opt, y=cue_strat)) +
geom_abline(slope=1, intercept=c(0), color='grey') +
labs(x='Optimum CUE (max B w/o compition)', y='Strategic CUE (resistant to invasion)')
#test.parm <- as.list(stratParm.df[1,])
#attach(test.parm)
#detach(test.parm)Explore drivers of strategic and optimum differences
diff.df <- subset(stratParm.df, validPools_opt & validPools_strat,
select=setdiff(names(stratParm.df), c('orderDiff_mh')))
diff.df <- melt(diff.df, measure.vars=names(parm.ls))
ggplot(diff.df) + geom_violin(aes(x=as.factor(value), y=(cue_opt-cue_strat)/cue_opt)) +
facet_wrap(~variable, scale='free')
We are particullary interested in differences in inputs
validParms <- subset(stratParm.df, validPools_strat & validPools_opt)
I.arr <- seq(0.1, 10, by=0.3)
Igrad <- ddply(unique(validParms[,c("b", "vmax", "kmin", "k_v_slope", "h", "m")]),
c("b", "vmax", "kmin", "k_v_slope", "h", "m"), function(xx){
ans1 <- data.frame()
for(I in I.arr){
opt_results <- optimum_cue(xx$b, xx$vmax, xx$kmin,
xx$k_v_slope, xx$h, xx$m, I)
names(opt_results) <- paste(names(opt_results), '_opt', sep='')
strat_results <- strategic_cueI(xx$b, xx$vmax, xx$kmin,
xx$k_v_slope, xx$m, xx$h, I)
ans1 <- rbind.fill(ans1,
as.data.frame(c(xx[c("b", "vmax", "kmin",
"k_v_slope", "h", "m")],
I=I, opt_results, strat_results)))
}
#cat('dim ans1:', dim(ans1), '\n')
ans2 <- data.frame()
for(I in I.arr*rnorm(length(I.arr), mean=1, sd=0.05)){
opt_results <- optimum_cue(xx$b, xx$vmax, xx$kmin, xx$k_v_slope,
xx$h, xx$m, I)
names(opt_results) <- paste(names(opt_results), '_opt', sep='')
strat_results <- strategic_cueI(xx$b, xx$vmax, xx$kmin, xx$k_v_slope,
xx$m, xx$h, I)
ans2 <- rbind.fill(ans2,
as.data.frame(c(xx[c("b", "vmax", "kmin",
"k_v_slope", "h", "m")],
I=I, opt_results, strat_results)))
}
#cat('dim ans2:', dim(ans1), '\n')
if(nrow(ans1) != nrow(ans2)){
print(xx[c("b", "vmax", "kmin", "k_v_slope", "h", "m"),])
return(NULL)
}
dCUEdI <- data.frame(dcue_opt = (ans1$cue_opt - ans2$cue_opt),
dcue_strat = (ans1$cue_strat-ans2$cue_strat),
dB_opt = ans1$B_opt-ans2$B_opt,
dC_opt = ans1$C_opt-ans2$C_opt,
dB_strat = ans1$B_strat - ans2$B_strat,
dC_strat = ans1$C_strat - ans2$C_strat,
dI = (ans1$I-ans2$I))
#cat('dim dCUEdI:', dim(dCUEdI), '\n')
return(cbind(ans1, dCUEdI))
})
Igrad$validParam <- with(Igrad, B_opt+C_opt > 10 &
B_opt+C_opt < 500 &
B_opt/(B_opt+C_opt)> 0.001 &
B_opt/(B_opt+C_opt) < 0.15 &
B_strat+C_strat > 10 &
B_strat+C_strat < 500 &
B_strat/(B_strat+C_strat)> 0.001 &
B_strat/(B_strat+C_strat) < 0.15)plot.df <- melt(Igrad, measure.vars=names(Igrad)[grepl('_(opt)|(strat)', names(Igrad))])
temp <- as.data.frame(matrix(unlist(strsplit(as.character(plot.df$variable), '_|\\.')), ncol=2, byrow=TRUE))
names(temp) <- c('variable', 'type')
plot.df$variable <- NULL
plot.df <- cbind(plot.df, temp)
plot.df$isGrad <- grepl('d', as.character(plot.df$variable))
ggplot(subset(plot.df, validParam & isGrad)) + geom_hex(aes(x=I, y=value/dI)) + facet_grid(variable~type, scales='free') + scale_y_log10()## Warning: Transformation introduced infinite values in continuous y-axis## Warning: Removed 70820 rows containing non-finite values (stat_binhex).
ggplot(subset(plot.df, validParam & !isGrad)) + geom_hex(aes(x=I, y=value)) + facet_grid(variable~type, scales='free')
Numerical validation code
Let’s pick some points where the gap between the Optimum CUE and Strategic CUE is small (\(<20\) percent) and large (\(>20\) percent).
bigGap.df <- subset(stratParm.df, validPools_opt & validPools_strat & abs(cue_strat-cue_opt) > 0.2)
smallGap.df <- subset(stratParm.df, validPools_opt & validPools_strat & abs(cue_strat-cue_opt) < 0.2)
numParm.df <- rbind.fill(bigGap.df[sample.int(nrow(bigGap.df), size=2),],
smallGap.df[sample.int(nrow(smallGap.df), size=2),])
numParm.df$index <- 1:nrow(numParm.df)
pander(numParm.df[,c(names(parm.ls), c('cue_opt','C_opt', 'B_opt', 'cue_strat','C_strat', 'B_strat'))])| b | vmax | kmin | k_v_slope | h | m | I | cue_opt | C_opt | B_opt |
|---|---|---|---|---|---|---|---|---|---|
| -1 | 10 | 100 | 100 | 0.1 | 1e-04 | 0.1 | 0.9394 | 56.45 | 0.8864 |
| 1 | 10 | 100 | 10 | 1 | 0.001 | 10 | 0.899 | 314.4 | 8.707 |
| -10 | 1 | 10 | 100 | 0.01 | 0.001 | 0.1 | 0.2828 | 23.72 | 2.158 |
| 0.1 | 10 | 100 | 1000 | 0.1 | 1e-04 | 1 | 0.9798 | 289.5 | 9.514 |
| cue_strat | C_strat | B_strat |
|---|---|---|
| 0.697 | 22.91 | 0.681 |
| 0.6162 | 71.39 | 6.118 |
| 0.1717 | 13.42 | 1.487 |
| 0.8687 | 134.8 | 8.57 |
Validate optimium CUE
opt.validate <- ddply(numParm.df, c('index', 'C_opt', 'B_opt'), function(xx){
cue <- seq(0, 1, length=100)
ans <- steadyState(xx$b, xx$vmax, xx$kmin, xx$k_v_slope, xx$h, xx$m, xx$I, cue)
return(data.frame(cue=cue, C=ans$C, B=ans$B))
})
plot.df <- melt(subset(opt.validate, B>0 & C>0), measure.vars=c('C', 'B'))
ggplot(plot.df) + geom_line(aes(x=cue, y=value, group=index)) +
geom_vline(data=numParm.df, aes(xintercept=cue_opt)) +
facet_grid(variable~index, scales='free')
Validate strategic CUE
compitition.model <- function(t, y, parms){
C <- y[1]; Bn <- y[2]; Bi <- y[3]
ans <- with(parms,
c(dC = I - m*C - C*(vn*Bn/(kn+C) + vi*Bi/(ki+C)),
dBn = cuen*vn*Bn*C/(kn+C) - hn*Bn,
dBi = cuei*vi*Bi*C/(ki+C) - hi*Bi)
)
return(list(ans))
}
runInvasions <- ddply(numParm.df, c('index'), function(xx){
parm <- as.list(xx[, c('b', 'vmax', 'kmin', 'k_v_slope', 'm', 'I')])
parm$hn <- xx$h
parm$hi <- xx$h
cueCombo <- expand.grid(cuen=seq(1, 100, length=50)/100,
rel_cuei=rnorm(10, sd=0.1))
cueCombo$cuei <- cueCombo$cuen+rnorm(nrow(cueCombo), sd=0.1)
cueCombo <- cueCombo[cueCombo$cuei > 0,]
cueCombo$vn <- with(parm, cue_v_tradeoff(b=b, vmax=vmax, cueCombo$cuen))
cueCombo$kn <- with(parm, v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=cueCombo$vn))
cueCombo$vi <- with(parm, cue_v_tradeoff(b=b, vmax=vmax, cueCombo$cuei))
cueCombo$ki <- with(parm, v_k_tradeoff(kmin=kmin, k_v_slope=k_v_slope, v=cueCombo$vi))
invade.df <- ddply(cueCombo, c('cuen', 'cuei'), function(cuePairs){
comboParm <- c(parm, cuePairs)
preInvade <- with(comboParm, steadyState(b, vmax, kmin, k_v_slope, h=hn, m, I, cue=cuen))
if(all(preInvade > 0)){
y0 <- list(C=preInvade$C, Bn=preInvade$B, Bi=preInvade$B*0.1)
invasion <- lsoda(y=unlist(y0), times=c(1, 7, 30, 365, 365*10), func=compitition.model, parms=comboParm)
names(y0) <- paste(names(y0), '0', sep='')
if(nrow(invasion) == 5){
finalPools <- as.data.frame(c(comboParm, y0, as.list(invasion[nrow(invasion),])))
}else{
finalPools <- as.data.frame(comboParm)
}
}else{
finalPools <- as.data.frame(comboParm)
}
})
return(invade.df)
})ggplot(runInvasions) + geom_point(aes(x=cuen, y=cuei, color=Bn > Bi)) +
geom_vline(data=numParm.df, aes(xintercept=cue_strat)) +
labs(x='Native CUE', y='Invader CUE', title='Winner of 10 year invasion') +
facet_wrap(~index)